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ABSTRACT 

In this work we address the problem of simultaneus multi-frequency detection of 
extragalactic point sources in maps of the Cosmic Microwave Background. We apply 
a new linear filtering technique, the so called 'matched matrix filters', that incorpo- 
rates full spatial information, including the cross-correlation among channels, without 
making any a priori assumption about the spectral behaviour of the sources. A sub- 
stantial reduction of the background is achieved thanks to the optimal combination 
of filtered maps. We describe in detail the new technique and we apply it to the de- 
tection/estimation of radio sources in realistic all-sky Planck simulations at 30, 44, 
70 and 100 GHz. Then we compare the results with the mono-frequencial approach 
based on the standard matched filter, in terms of reliablity, completeness and flux ac- 
curacy of the resulting point source catalogs. The new filters outperform the standard 
matched filters for all these indexes at 30, 44 and 70 GHz, whereas at 100 GHz both 
kind of filters have a similar performance. We find a noticeable increment of the num- 
ber of true detections for a fixed reliability level. In particular, for a 95% reliability 
we practically double the number of detections at 30, 44 and 70 GHz. 

Key words: methods: data analysis - techniques: image processing - radio contin- 
uum: galaxies - cosmic microwave background - surveys 



1 INTRODUCTION 

From the search of extrasolar planets to the study of ac- 
tive galactic nuclei, one of the most common task in all the 
branches of Astronomy is the detection of faint pointlike ob- 
jects. Such objects have angular sizes that are smaller than 
the angular resolution of the telescopes that are used to ob- 
serve them, and therefore they are usually referred to as 
point sources. 

A case of particular interest is the detection of extra- 
galactic point sources (EPS) in maps of the Cosmic Mi- 
crowave Background (CMB). EPS are known to be a rel- 
evant source of contamination for CMB studies, specially 
at small angular scales, where they hamper the estima- 
tion of CMB angular power spectrum both i n temperature 



necessary to detect and remove as many extragalactic point 
souces as possible. 

Moreover, in the frequency range spanned by CMB 
experiments the properties of EPS are poorly studied. 
Only very recently the Wilkinson Micro wave Anisotropy 
Probe (WMAP) satellite ([Bennett et all l2Q03h has per- 
mitted the obtention of the first all-sky complete point 
source catalogs above ^ 0.8-1 Jy in the 23-94 GHz 
of frequencies (|Bennett et al. 20031: iHinshaw et al.1 



tion ol UM-b angular power spectrum both m tem perature 
dToffolatti et al.1 Il998l: Ide Zotti et all Il999l : iHobson et al 



1 9991: Ide Zotti et alJl2005h and in polarization ([Tucci et al 



12004 l2005h . Therefore, for the sake of CMB analysis it is 
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l2007t iLopez-Caniego et al.1 120071 : IChen fc Wright! 120071: 
[ Wright et al.ll2QQ8h . The upcoming Planck mission (iTauberl 
I2005T ) will allow us to extend these catalogs down to 
lower flux limits and up to 857 GHz. There is interest- 
ing physics to be probed in this frequency range. The new 
EPS catalogs provided by next generation CMB experi- 
ments will not only allow us to follow the behaviour of 
source counts from existing catalogs to microwave frequen- 
cies, but also to study source variability and to discover 
rare objects such as inverted spectrum radio sources, ex- 
treme gigahertz peaked spectrum (GPS) sources and high- 
redshift dusty galaxies (see for example the Planck Blue- 



2 Herranz et al 



book (Th e Planck Collaboration! l2006h for a brief, yet com- 
prehensive review of the rich phenomenology of EPS at 
microwave frequencies). Thus, the task of detecting point 
sources is important not only from the point of view of CMB 
science but also from the point of view of extragalactic As- 
tronomy as well. 

Let us consider a single image taken at a given 
wavelength. Then the problem consists on how to detect a 
number of objects, all of them with a common waveform 
that is generally considered to be well known (basically, 
the shape of point sources is that of the beam) but with 
unknown positions and intensities, that are embedded 
in additive noise (not necessarily white). In the field of 
micro w ave A st ronomy, wavelet tech ni ques ([Vielva et al 
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2005) have proved to be useful. The 



common feature of all these techniques is that they rely 
on the prior knowledge that the sources have a distinctive 
spatial behaviour (i.e. a known spatial profile, plus the fact 
that they appear as compact objects as opposed to 'diffuse' 
random fields) that helps to distinguish them from the 
noise. Detection can be further improved by including prior 
information about the sources, i.e. some knowledge about 
their intensit y distribution, in th e fram e of a Bayesian 
form alism ([Hobson &; McLachla 3 120031 : ICarvalho et aU 

Most of the current and planned CMB experiments are 
able to observe the sky at several wavelengths simultane- 
ously. Multi-wavelength information makes it possible to 
separate different astrophysical components (as for example 
CMB from Galactic synchrotron emission) that have differ- 
ent spectral behaviour. Although mult i wavelength compo- 
nent separation techniques have been very succesful in sep- 
arating diffuse components (for a recent comparative review 
of several methods applied to sky sim ulations very simil ar to 
the ones we will use in this paper, see lLeach et al. 2008|), the 
detection of EPS has been usually attempted on a channel 
by channel basis. The reason for this is that EPS form a very 
heterogeneous population, constituted by a large number of 
objects with very different physical properties, and therefore 
it is impossible to define a common spectral behaviour for 
all of them. The whole situation remains somewhat unsatis- 
factory: on the one hand, the channel- by- channel approach 
based on the spatial behaviour works fine, but a valuable 
fraction of the information that multi-wavelenght experi- 
ments can offer is wasted this way. On the other hand, stan- 
dard component separation techniques based on spectral di- 
versity have problems when dealing with the very heteroge- 
neous EPS components. 

Thus, multi-wavelenght detection of EPS in CMB im- 
ages remains a largely unexplored field. In recent years, some 
attempts have been done in this direction. For example, 
iNaselskv et~aH (|2002h combined simulated multi-wavelenght 
maps ir order to increase the avera ge signal to noise ratio 
of point sources. In a similar way, IChen fe Wrightl (|2007T ) 
and I Wright et all (|2008T ) use combinations of the WMAP W 
and V bands in order to produce a CMB-free map in which 



to better detect the elusive radio galaxies. Note that, in any 
case, combined 'clean' maps are suitable for detecting more 
sources but not for performing accurate photometry, unless 
the spectral index of all the sources is known in advance. 

An intermediate approach is to design filters that are 
able to find compact sources thanks to their distinctive 
spatial behaviour while at the same time do incorporate 
some mult i wavelength information, without pretending to 
achieve a full component separation and without assuming 
a specific spectral behaviour for the sources. Very recently, 
the authors have proposed a new techniq ue, based on the 
so-ca lled 'matched matrix filters' fMTXF, iHerranz &; Sanzl 
I2008T ) , that goes in this direction. The basic underlying ideas 
of the new method are: 

• When a source is found in one channel, it will be also 
present in the same position in all the other channels. 

• The spatial profile of the sources may differ from chan- 
nel to channel, but it is a priori known. 

• The second order statistics of the background in which 
the sources are embedded is well known or it can be directly 
estimated from the data by assuming that point sources 
are sparse. This knowledge about the second order statistics 
(namely, the background's power spectrum for each channel 
and the its correlations among the different channels) will 
be used to increase the signal to noise ratio of the sources. 

• We want to perform accurate photometry of the sources 
at every one of the frequencies covered by the experiment, in- 
dependently of what is the spectral behaviour of any source 
in particular. 

In lHerranz~fc Sand (|2008h the authors presented the new 
methodology and demonstrated its potential utility with a 
few toy simulations. In this paper, we will study its appli- 
cability to real CMB experiments by applying it to realistic 
simulations of the whole sky as will be observed by the up- 
coming Planck mission. We will focus on the particular case 
of the detection of radio sources in the four lower frequency 
Planck channels (33-100 GHz), comparing the performance 
of the new filters with the performance of the well-stablished 
standard matched filters. In section [2] we will summarize the 
foundations of the theoretical formulation of matched ma- 
trix filters. In section [3] we will describe the Planck simula- 
tions that we use to test the method and we will outline the 
main features of the code we have developed for its imple- 
mentation. The results of the exercise will be commented in 
section 2] Finally, in section [5] we will draw some conclusions. 



2 MATCHED MATRIX FILTERS 

The deriv ation of the ma t ched matrix filters is fully de- 
scribed in Herranz & Sana (20080 However, for reasons of 
clarity, we will reproduce here the main ideas of that paper 
here. 



2.1 Data model 

Let us consider a set of N two-dimensional images (chan- 
nels) in which there is an unknown number of point sources 



1 For economy, in the following we will occasionally refer to the 
matched matrix niters just as 'matrix niters' 
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embedded in a mixture of instrumental noise and other as- 
trophysical components. Without loss of generality, let us 
consider the case of a single point source located at the ori- 
gin of the coordinates. Our data model is 



D k (x) = s k (x) + n k (x) , 



(1) 



where the subscript k = 1 , . . . , N denotes the index of the 
image. The term Sk(x) denotes the point source, 

s k (x) = A k r k (x) , (2) 

where Ak is the unknown amplitude of the source in the k th 
channel and Tk (x) is the spatial profile of the source (which is 
assumed to be known) and satisfies the condition Tfc(O) = 1. 
The term rik(x) in equation (pQ) is the generalized noise in 
the k th channel, containing not only instrumental noise, but 
also CMB and all the other astrophysical components apart 
from the point sources. Let us suppose the noise term can 
be characterized by its cross-power spectrum: 



(rife (q)ni (q')) = P k i (q) S 2 (q - q' 



(3) 



where P = (Pki) is the cross-power spectrum matrix and 
the symbol * denotes complex conjugation. From now on, 
we assume that the noise has zero mean: 



(rife (a?)) - 0. 



2.2 Filtering with matrices of filters 



(4) 



Since we are interested in doing accurate photometry in each 
one of the N available channels, we are bound to produce 
N different processed maps. Therefore, we are looking for a 
transformation that starts with N input channels and ends 
with other N processed maps where a) point sources are 
easier to detect and b) the amplitudes Ak are preserved. Be- 
sides, since we intend to use some mult i- wavelength informa- 
tion, we are interested in making all the N input channels 
intervene in the elaboration of the any one of the output 
maps. One possibility is to define a set of N x N filters 
^ki(x) such that the N combined quantities 



w k {x) 
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i 
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Di (i 



*fc( (q) D, (q) 



(5) 



are our processed maps. The last term of the equation is just 
the expression of the filterings in Fourier space, being q the 
Fourier mode and ^ik(q) and Di(q) the Fourier transforms 
of ^ki(x) and Di(x), respectively. 

We intend to use the combined filtered image Wk(x) 
as an estimator of the source amplitudes Ak for all k = 
1,...,N. Thus, the filters ^ki must satisfy the condition 
that the k th filtered, combined image at the position of the 
source is, on average over many realizations, an unbiased 
estimator of the amplitude of the k th source. The other re- 
quirement we will ask the processed maps is that the signal 
to noise ratio of the sources is increased with respect to the 
input maps. In other words, we want estimators Wk to be 
not only unbiased, but efficient as well. Therefore, we need 
to minimize the variance a Wk of the combined filtered image. 

The set of filters that minimize the variance cr Wk for all 
k while keeping the individual amplitudes Ak constant for 




Figure 1. Mask used for the analysis. The mask covers 14.91% 
of the sky, including the Galactic Plane plus some other densely 
contamined areas of the sky such as the Magellanic Clouds, the 
Ophiuchus Complex and the Orion/Eridanus Bubble. 



all the point sources, independently of their frequency de- 
pendence, can be shown to be given by the matrix equation: 



ijr* = FP" 1 , 

where 

F = (F kl ), P = (P k i), A = (Afcj), H = (H kl ), 
and where 



(6) 



Fki 
A 



Hkl = J dqr k (q) P^niq) . 



(8) 



The set of filters we have developed naturally assume a 
structure that is best expressed in the form of a matrix equa- 
tion, hence the denomination of 'matched matrix filters' 



2.3 Properties of matched matrix filters and some 
particular cases 

2.3.1 A single image 

If the number of images is N — 1 it is easy to show that the 
matrix of filters contains a single element, which is the com- 
plex conjugate of the standard matched filter. For circularly 
symmetric source profiles, the filter is real- valued and the 
resulting filter is exactly the same as the standard matched 
filter. 



2.3.2 Uncorrelated noise 

From equations (|6} and J5J) it is straightforward to show that 
for the particular case where the noise is totally uncorrelated 
among channels the matrix of filters defaults to a diagonal 
matrix whose non-zero elements are the complex conjugates 
of the standard matched filters that correspond to each input 
channel. When the source profiles are circularly symmetric 
the filters are real- valued and the whole process is equivalent 
to filter each channel independiently with the appropiate 
matched filter. 
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Figure 2. Standard matched filters, in Fourier space, for the four 
channels in one of the patches considered. 



2.3.3 The 2 x 2 case 

In this case, the explicit form of the MTXF is 
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where z, j = 1, 2 and 0j are closely related to the standard 
matched filters ? MF : 



£.3.^ Relative gain of the different channels 

An useful quantity that can be used to assess the perfor- 
mance of a filter is its gain factor, that is, the increment of 
the signal to noise ratio that a fiducial source experiments 
thanks to the application of a filter. Of course, gain is not a 
universal quantity of the filter, since it depends on the sta- 
tistical properties of the background noise, but given a par- 
ticular image the relative gain of two filters provides a good 
intuitive measurement of their respective performances. For 
unbiased filters (i.e. filters that preserve the amplitude of the 
signal) the gain is just the ratio of noise dispersions before 
and after filtering: 



G 



a 

0~ w 



(13) 



where a is the noise dispersion before filtering and a w is the 
noise dispersion after filtering. 

When comparing the performances of matched 
matr ix filters with thos e of standard matched fil- 
ters, lHerranz"^T Sand (2008) observed an apparent trend: in 
all the few cases they studied, the gain factor of the matrix 
niters was similar to the one of matched filters for at least 
one of the channels, and significantly higher for the others. 
Moreover, the channel with the lower gain was the one with 
the worse signal to noise ratio after filtering. At the mo- 
ment, it was not clear if this was a universal trend of just 
a result obtained by chance given the low number of cases 
under studied in that work. As we will see below, the results 
of this work seem to support this observed trend. 

There is a qualitative argument that sheds light over 
this phenomenon: for the sake of simplicity let us consider 
two channels and identical profiles n = T2. If one assumes 
that Pn < P22 then the equations (|10lll|) lead to 



< cr u 



MF MF 



(14) 



for the MTXF and standard MF, respectively. Therefore, the 
original map with less variance gains more with either filter. 
On the other hand, assuming that P12 <C P11 < P22, taking 
into account equations (|10llip and the Schwartz inequality, 
one obtains 



^ CT 7( 



i.e. the MTXF outperform the standard MF. 



(15) 
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We remark that the MTXF is a non-symmetric matrix in 
general. The variances of the two filtered maps are 



iV6n' 
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From the last equation, one obtains 
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3 APPLICATION TO PLANCK RADIO 
SOURCES 

In this section, we will describe an application of the new 
mult i wavelength filtering technique to realistic simulations 
of the sky as it will be observed by the Planck mission. As 
a example, we will focus on the blind detection of extra- 
galactic radio sources. Therefore we will simulate the 30, 
44, 70 and 100 GHz Planck channels. Note that although 
a simultaneous mult i wavelength filtering of the nine Planck 
channels is relatively easy to do (it would imply the use of 
a 9 x 9 matrix of niters, that can be easily handled by any 
desktop computer), the addition of infrared channels would 
probably help little to the detection radio sources. A intu- 
itive reason for this is that the sources that are bright in the 
infrared range are not necesarily the same that dominate 
the radio range, and vice versa. Therefore, for this work 



Multi-frequency detection of point sources in CMB maps 5 




Figure 3. The different elements of the 4x4 matrix of niters, in Fourier space, for the same patches as in figure [2] 



we prefer to treat the radio sources separately. Besides, the 
low- dimensional 4x4 application we describe here is more 
adequate to illustrate the technique, allowing us to show 
the details of the method with a relatively small number of 
plots. 



3.1 The simulations 

Sky simulations are bas ed on the Planck Sky ModeQ (PSM, 
iDelabrouille et aDl2008l in preparation), a flexible software 
package developed by Planck WG2 for making predictions, 
simulations and constrained realisations of the microwave 
sky. 

The CMB sky is based on a Gaussian realisation assum- 
ing the WMAP best-fit Ce at higher multipoles. 

The Galactic emission is described by a three compo- 
nent model of the interstellar medium comprising free- free, 
synchrotron and dust emissions. Free-free emission is based 



2 http://www.apc.univ-paris7.fr/APC_CS/Recherche/ 
Adamis/PSM / psky-en.php 



on the model of [Dickinson et aL I (|2003h assuming an elec- 
tronic temperature of 7000 K. The spatial structure of the 
emission is estimated using a Ha template corrected for dust 
extinction. Synchrotron e mission is based on a n extrapola- 
tion of the 408 MHz map of Has lam et al.l (|1982T ) from which 
an estimate of the free- free emission was removed. A limi- 
tation of this approach is that this synchrotron model also 
contains any dust anomalous emission seen by WMAP at 
23 GHz. The thermal e mission from inter s tellar dust is esti- 
mated using model 7 of iFinkbeiner et al.1 (|l999h . 

Point sources are modeled with two main categories: 
radio and infra-red. Simulated radio sources are based on 
the NVSS or SUMSS and GB6 or PMN catalogues. Mea- 
sured fluxes at 1 and/or 4.85 GHz are extrapolated to 
PLANCK frequencies assuming a distribution in fiat and 
steep populations. For each of these two populations, the 
spectral index is randomly drawn within a set of values 
compatible with the typical average and dispersion. Infrared 
sources are based on the IRAS catalogue, and modelled as 
dusty galaxies. In addition, the emission of a large num- 
ber of blended infrared galaxies, not present individually 
in the IRAS catalogue, is simulated to model the Far In- 
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frared Background (| Gonzalez- Nuevo et al.l l2005h . We also 
include in the model a map of thermal SZ spectral distortion 
from galaxy clusters, based on a cluster catalogue randomly 
drawn using a mass-function compatible with present-day 
observations. 

3.2 The code 

In this work we have a used a code that reads in four all- sky 
FITS maps with HEALPix (jGorski et all l2Q05h resolution 
parameter NSIDE=1024, one per frequency between 30 and 
100 GHz. Second, according to some parameters given in 
an input file, the code uses the CPACK libraries to divide 
the sky in a sufficient number of overlapping flat patches 
such that the 100% of the sky is covered. In this pixelization 
scheme, we have produced 371 patches per frequency, 14.656 
square degrees and 256 x 256 pixels each. Then the code 
proceeds to run either the MF or the MTXF algorithms on 
every set of four patches corresponding to the same region 
in the sky. Afterwards, once the optimally filtered image has 
been produced, the code looks for maxima in it producing a 
subcatalogue of detections. Finally, a combined catalogue is 
produced, removing possible repetitions inside a 1 FWHM 
radius. 



4 RESULTS 

As described above, we have divided t he sky into square 
flat patches by projecting the HEALPix ([Gorski et al.| [2005) 
maps into the tangent plane at a set of coordinates that are 
regularly distributed on the sphere. For each resulting patch 
we have four square images, one per frequency channel (30, 
44, 70 and 100 GHz). Then we have filtered simultaneously 
the four images with the filters J6]) specifically calculated for 
that region of the sky. In paralel, we have filtered each one of 
the images separately with its corresponding matched filter. 
Therefore, for each input image we have two output filtered 
images, one obtained with the standard matched filter and 
other obtained by the matrix filters as in equation J5}. Then 
we have applied the same thresholding detection criterion to 
the two cases: for any given filtered map, we have selected 
all the peaks that have at least three connected pixels with 
flux above a given number of times the a levejj of the fil- 
tered map. Unless otherwise noted, all the plots that will 
be shown in this section will refer to detections above the 
5cr detection threshold. The detections thus obtained for the 
different patches have been combined into a single whole-sky 
catalogs for each of the four frequency channels and the two 
filtering schemes. 

Figure [T] shows the Galactic mask we apply for the anal- 
ysis of the results. The mask is similar to the WMAP Kp2 
mask and it covers a highly contamined region around the 
Galactic plane plus a set of irregular areas that mask other 
highly contamied areas of the sky such as the Magellanic 
Clouds, the Ophiuchus Complex and the Orion/Eridanus 

3 http:/ /astro. ic.ac.uk/~mortlock/cpack/ 

4 Where a is the standard deviation of the filtered map, excluding 
the region of the borders of the image; note that this a level 
corresponds to a different flux threshold for different filters and 
for different regions of the sky 



Bubble. In total, we are masking 14.91% of the pixels of the 
sky. 

4.1 Detail of a single patch 

Before discussing the results for all the sky outside the mask 
just described above, let us illustrate the qualitative func- 
tioning of the filters taking as a example just one sky patch. 
In the first row of plots of figure |4] we show the aspect of 
sky in the first of the patches we have studied, centered in 
the Galactic North Pole. It is a region of the sky with a very 
low contamination from Galactic emission, with two point 
sources that are clearly visible to the naked eye (at least at 
30 GHz) , but many others are hidden amid the diffuse com- 
ponents. EPS alone are shown in the second row of plots of 
the figure. 

Figure [2] shows the matched filters, in Fourier space, 
for the four different channels shown in figure 2) Figure [3] 
shows the 16 elements of the corresponding matrix of filters, 
also in Fourier space. Note that the filters in the diagonal 
look roughly similar to the matched filters whereas the off- 
diagonal elements are quite different. This can be intuitively 
explained in the following way: the diagonal element ^kk is 
designed to produced a maximum contribution of the source 
profile Tfc in the map Dj~ whereas the ^ki, I ^ k element is 
designed to produce a minimum contribution of the source 
profile ti in the map Dk . This way the off-diagonal elements 
of the filtering contribute to reduce noise but do not intro- 
duce bias in the determination of the fluxes in the k th map. 

The third and fourth rows of figure [4] show the output 
filtered patches for the matched filter and the matrix fil- 
ters, respectively. Note that for the 44 and 70 GHz channels 
the output matrix filtered maps look far cleaner than their 
matched filtered equivalents. For the 30 GHz channel the 
distinction is not so clear (the matrix filtered image looks 
cleaner, but some of the sources that are easily visible in the 
matched filtered image are apparently missing; we will see 
later that this is only a visual effect). Finally, for the 100 
GHz channel both filtered images look practically identical. 

The gain factors obtained for these images with the 
MTXF are [2.9,3.8,3.5,2.8] for the [30,44,70,100] GHz 
channels. The gain ratio between the MTXF and the 
MF are Gmtxf/Gmf = [1.38,1.52,1.49,1.00] for the 
[30, 44, 70, 100] GHz channels. Thus, the MTXF-filtered im- 
age with lower gain (100 GHz) is the one that is more simi- 
lar to the correspondent MF-filtered image. Besides, the 100 
GHz map is the one with higher variance before filtering. We 
will see in the next section that the all-sky results confirm 
this rule. 



4.2 All-sky results 

We can use the knowledge on the input EPS we have simu- 
lated to control the quality of our catalogs in terms of num- 
ber of true and false detections. The criterion we use to 
decide whether a detection is a true one or a false one is 
purely positional: an object of the catalog is considered a 
true detection if it is closer than a certain matching radius r 
of a input source and a false (spurious) detection otherwise. 
For this work we use a matching radius r = 2Ro , where Ro is 
the width of the Gaussian beam correspondig to the channel 
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Figure 4. One of the regions of the sky. The patch is centered in the North Galactic Pole. The first row of images, starting from the top, 
shows the input patches at 30, 44, 70 and 100 GHz (from left to right of the figure). The second row of images shows separately the EPS 
contribution to the input maps, for the same frequencies. The third row shows the patches after having been filtered with the standard 
matched filter corresponding to each frequency. The bottom row shows the combined filtered maps resulting from the application of the 
matrix filters. All maps are expressed in MJy/sr units. 



under study. This radius is smaller than one FWHM, so the 
criterion we use is quite stringent. Considerations about the 
flux matching will be made later in section 14.2.51 

4-2.1 Number of detections 

Figure [5] shows the number of detected sources in our 5cr 
catalogs that have true fluxes S > So, as a function of the 
flux threshold So. The number of detections obtained with 
the MF is shown by a solid black line, whereas the num- 
ber of detections obtained with the MTXF is shown by a 
dot-dashed blue line (colors are available only in the online 
version of the paper). 

Two different cases can be observed: the 100 GHz chan- 
nel and the other three channels. At the 100 GHz channel 
the performance of the two kinds of filters is almost identical. 
From 30 to 70 GHz and high and intermediate fluxes (^ 0.6 
Jy) the two methods detect very similarly. The MTXF curve 



runs slightly below the MF curve: the difference consists of 
a few high flux sources that are somehow missed by the 
MTXF. We will discuss this problem in more detail below. 

In the low flux region of the plots the number of de- 
tected sources stop growing and the curves reache a plateau. 
The knee point of the curves roughly indicates the detection 
limit of the 5cr catalogs. At 30, 44 and 70 GHz this turnover 
point occurs around 600 (400) mJy for the matched (matrix) 
filters. In other words, the MF reaches its detection limit at 
fluxes higher than the MTXF. Therefore the MTXF are able 
to go deeper and to detect many more faint sources than the 
MF. At 100 GHz both filters have their turnout points lo- 
cated around 400 mJy. 

We have checked that the number of 5cr detections ob- 
tained with the MF roughly agrees (taking into account the 
differen t sky coverage) with the results obtained in pr evious 
works (jLoDez-Caniego et all 120061 : iLeach et all l2008h that 
have made use of similar Planck simulations. 
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Figure 5. Number of true detections with true fluxes above a 
given flux value. Black solid line: detections with the matched 
filter. Blue dot-dashed line: detections with the matrix filters. 
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Figure 6. Completeness diagram. Black solid line: matched fil- 
ter. Blue dashed line: matrix filters. The 100% and the 95% com- 
pleteness levels are marked by a red dotted and a green dotted 
horizontal lines, respectively. 



4-2.2 Completeness 

Figure [6] shows the completeness level as a function of the 
flux of the 5cr catalogs obtained with the MF and the MTXF. 
Here the completeness is defined, as usual, as the ratio be- 
tween the number of recovered true sources and the total 
number of input sources over a given flux limit. The 95% 
completeness level is marked by a horizontal green dotted 
line. The 95% completeness fluxes are, for the matched (ma- 
trix) filters, the following: at 30 GHz, 610 (540) mJy; at 44 
GHz, 460 (340) mJy; at 70 GHz, 390 (270) mJy and, at 100 
GHz, 270 (270) mJy. Three issues about figure [6] deserve 
detailed comments: 

If compared to the data in Table 2 
of lLopez-Caniego et al.l ([2006), our present complete- 



ness limits are higher than theirs. This is due to two causes: 
on the one hand, in this work we include regions of the 
sky that are closer to the Galactic plane and therefore 
more contamined (and therefore, the 5cr detection threshold 
corresponds to a higher flux), where it is more likely to 
loss some sources, and, on the other hand, our present 
detection/selection criterion requires to have peaks with at 
least three connected pixels. This more stringent criterion 
helps to reduce the number of spurious detections, but 
ocasionally makes us loss some true sources as well. 

Regarding the different channels, for the case of 30, 44 
and 70 GHz the MTXF show a better completeness level 
at intermediate and low fluxes. This is because MTXF do 
amplify better the point source signal with respect to the 
foregrounds and therefore they can reach lower detection 
limits. At 100 GHz, both kinds of filtering lead to the same 
completeness levels. 

Both methods do miss some bright sources, even at 
fluxes > 1 Jy. For example, at 30 GHz both the MF and 
the MTXF fail to detect 4 sources with fluxes > 1 Jy. Three 
of them are common for the two filters and in all the three 
cases are sources that are in heavily contamined regions at 
low Galactic latitude, very close to the border of the mask. 
Apart from these three missing sources, the MF misses a 
source that is detected by the MTXF and the MTXF miss a 
source that is detected by the MF. The first one corresponds 
to a source close to the LMC and a few pixels away from 
our mask. The second case corresponds to a source that is 
2.69 Jy source that is only 35 arcmin away from a detected 
3.23 Jy source. Although the distance between both sources 
is slightly larger than the matching radius we used for this 
channel, the MTXF is unable in this case to resolve the two 
sources individually. 

This suggest that the mask we have used is good enough 
to avoid most spurious detections, but it may be insufficient 
to guarantee completeness. Besides, a problem of blending 
with the MTXF arises when two high flux sources lie very 
close one to another. Although this situation is not frequent, 
this may explain the slight loss of performance of the MTXF 
in figure [5] for high and intermediate fluxes. 



4.2.3 Reliability 

Another interesting indicator of the performance of the fil- 
ters is the reliability of the catalogs obtained with them. Let 
it be Nd(y) the number of true detections above a certain 
detection threshold v and N s (y) the number of spurious de- 
tections (false alarms) above the same threshold. Then we 
define the reliability in the usual way 



N d {y) 



(16) 



N d (v) + N s (v)' 

The reliability of our 5cr catalogs is well above the 95% 
level for the two filtering schemes and the four channels we 
considered. In order to show the performance of the filter 
for lower reliability levels, in figure [71 we have gone deeper 
in the detection, down to the 3a level. The red dotted line 
in the figure [71 shows the r = 95% reliability level. For the 
30, 44 and 70 GHz channels the MTXF allow us to go to 
lower detection thresholds than the MF for a fixed required 
reliability level. For the 100 GHz channel the situation is 
the opposite, but the difference is small. Table [1] shows the 
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Figure 7. Reliability of the catalogs as a function of the limiting 
a threshold for the two filters and the four channels. Black solid 
line: MF. Blue dot-dashed line: MTXF. The red dotted line shows 
the 95% reliability level. 



threshold limits at the 95% reliability level and the corre- 
sponding number of true detections for the two niters and 
four channels considered. The improvement of the MTXF 
with respect to the MF is significant for the 30, 44 and 70 
GHz channels. 



4.2.4 Receiver operating characteristics 

Yet another way to comparatively study the performance of 
two detectors is the so-called receiver operating characteris- 
tic (ROC), or simply ROC curve. ROC curves are profusely 
used in detection theory because they provide a direct and 
natural way to relate the costs/benefits of the decision mak- 
ing associated to the detection process. Let us consider the 
two following quantities: the true positives ratio (TPR) is 
denned as 



TPR: 



N d 
NT' 



(17) 



where, as before, Nd is the number of true detections (true 
positives) obtained for a certain detection threshold and NT 
is the total number of objects (in this case, simulated point 
sources) in the data set. The TPR is related, but not equiv- 
alent, to the completeness denned above (the number in the 
denominator of the completeness depends as well on the de- 
tection threshold, but NT does not). 

The other quantity of interest is the spurious detection, 
false alarm or false positive ratio (FPR) : 



FPR: 



Ns_ 
FT' 



(18) 



where N s is denned as in section 14.2.31 and FT is the total 
number of candidates (in our case, peaks in the filtered im- 
ages) that can be identified or not as 'detections' (true or 
false) by the detector. Therefore, the is an indirect relation 
between FPR and the reliability, but note that the quan- 
tity in the denominator in eq. (|18[) does not depend on the 
threshold. Both quantities, TPR and FPR, take values in 



Figure 8. ROC curves for the filtering schemes and the four 
channels considered. Black solid line: matched filter. Blue dashed 
line: matrix filters. 



the interval [0, 1] and they are called the operating charac- 
teristics of the detector. The TPR can be directly associated 
to the power of the detector and the FPR is related to its 
significance. 

ROC curves are constructed by plotting the fraction of 
true positives (TPR) vs. the fraction of false alarms (FPR). 
They convey at a single glance the same information that 
can be found by taking together the figure plus a set of 
tables akin to table [1] obtained at different reliability levels. 
For any fixed false alarm ratio the ROC curve tells the (nor- 
malized) number of true detections we have. ROC curves 
facilitate the comparison between two or more detectors (in 
our case, niters): the curve that lies above in the plot is 
closer to the optimal performance than the curves below. 

Figure [8] shows the ROC curves for the different cases 
and channels under consideration. The black solid line cor- 
responds to the matched filter and the blue dashed line cor- 
responds to the matrix niters (colors are available in the 
online version of the paper and on request to the authors). 
The number NT used for the plot is the total number of 
sources that are outside the masked area of the sky and 
that have fluxes above 150 mJy. The number FT is the to- 
tal number of maxima found in the filtered images above 
the mimimum considered threshold. The plots are made for 
a detection thresholds in the interval a £ [3, 10]. For this 
range of detection thresholds, the ROC curve corresponding 
to the matrix niters is clearly above the one corresponding 
to the matched niters for the 30, 44 and 70 GHz channels, 
meaning that if we fix any required spurious detection ratio 
we always have more true detections with the matrix niters 
than with the matched filters. For the 100 GHz channels the 
two curves run practically in parallel, but the matched fil- 
ters are slightly above the matrix filters. This is related to 
the missing sources problem described in section 14.2.21 



4-2.5 Flux estimation 

Figure [9] shows how both kinds of niters do recover the fluxes 
of the sources in all the considered cases. Circles represent 
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Table 1. Threshold c 95% limit required for a 95% reliability, and corresponding number N d ^ 95 % of detections for such a threshold. 
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Figure 9. Estimation of the flux for the two filtering techniques. 
Red circles: matched filter estimation. Blue crosses: matrix filter 
estimation. 
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Figure 10. Relative flux estimation error and its dispersion. 
Black lines: matched filters. Blue solid line plus dots: matrix fil- 
ters. For clarity, the lines corresponding to the matrix filter have 
been slightly displaced towards the right. 



fluxes recovered with the MF and crosses represent fluxes 
recovered with the MTXF. There is an excellent agreement 
between true input and recovered fluxes for all the cases and 
between fluxes obtained with the MF and the MTXF. At 
low fluxes, the MF estimates show the well-know selection 
Eddington bias before than the MTXF. This is particularly 
evident at 44 and 70 GHz, but also visible at 30 GHz. As 
seen before, the MTXF filtered maps have less noise than the 
maps filtered with the MF and therefore the selection bias 
appears at a lower level. This also manifests in the smaller 
dispersion of blue crosses around the true value with respect 
to the red circles. However, for very high fluxes it seems that 
MTXF tends to underestimate the flux of the sources. At 100 
GHz both filters lead to virtually the same results. 

Figure [10] shows the relative flux error for the two filters 
and the four considered channels. We define the relative flux 
error as 

e reZ = 100x (19) 

where S is the flux we estimate with the filters and So is 
the input flux. We take five flux bins between the flux of 
the faintest detected source and 2 Jy and for each bin we 
compute the average value and the dispersion of the relative 
error. The lowest bin is dominated by Eddington bias. This 
bias can be quite large, as is is the case for the MF at 44 
GHz (115% bias). The MTXF suffer less Eddinton bias than 
the MF (except, as usual, for the 100 GHz case, where the 
two filters behave very similarly). As mentioned before, the 



dispersion around the mean value is smaller for the MTXF. 
Finally, it is worth noticing the small negative bias suffered 
by both filters at high fluxes. This bias is not evident at 30 
GHz, but it is easier to detect at higher frequencies. MTXF 
seem to be more biased (—6%, —6% and —8% at 44, 70 and 
100 GHz, versus the MF biases that are -1%, -3% and 
—8% for the same frequencies). This kind of effect at high 
fluxes has been noticed be fore in related filtering applica- 
tions (Herranz et "all l2002h and it is usually attributed to 
the non-ideality of the pixelized data, that makes imperfect 
the process of calculation of the normalization of the filters. 
In any case, the bias is relatively small and can be calibrated 
by means of simulations. 



5 CONCLUSIONS 

Although there is a relatively large number of works in the 
literature devoted to the detection of point sources in CMB 
images, there is practically no one that addresses the prob- 
lem from the multi-frequency point of view. The reason is 
that each individual extragalactic point source has its own 
(a priori unknown) spectral behaviour and this makes it very 
hard to accomodate them in classic component separation 
schemes. 

In this work we apply a novel linear fil tering technique, 
the 'm atched matrix filters \ introduced by iHerranz &; Sanzl 
(2008). Matched matrix filters incorporate full spatial in- 
formation, including the cross-correlation among channels, 
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without making any a priori assumption about the spectral 
behaviour of the sources. The basic underlying idea is that 
in all the considered channels the sources do appear in the 
same unknown positions (but with different, unknown inten- 
sities) and have known spatial profiles given by the experi- 
ment point spread function, while some components of the 
background (i.e. CMB and Galactic emission) are correlated 
among channels. Then a clever linear filtering/combination 
of images can lead to a substantial reduction of the back- 
ground, and therefore to lower source detection thresholds. 
The resulting expresion of the filters takes form of a matrix 
given in terms of the cross-power spectra and the profile of 
the sources in the different maps. 

We describe in detail the formalism of the matched ma- 
trix filters, looking in detail into some particular cases of in- 
terest, and we apply them to the detection of radio sources 
in realistic all-sky Planck simulations at 30, 44, 70 and 100 
GHz. We use state-of-the-art simulations developed by the 
Planck WG2 Collaboration that include all the astrophysical 
components and instrumental specifications of the Planck 
mission. In order to tackle the strong non stationarity of the 
sky emission, we divide the sky into 371 14.656 x 14.656 over- 
lapping square degrees flat patches where we apply the two 
filtering techniques. For each patch and frequency we obtain 
a sub- catalog of detections. Then all the sub-catalogs corre- 
sponding to the same frequency are combined into an all-sky 
catalog. In order to compare with a well stablished mono- 
frequencial approach, we repeate the same process using the 
standard matched filter. 

We compare both methods in terms of reliablity, com- 
pleteness, receiver operating characteristics and flux accu- 
racy. We find that for the three lower frequency channels 
(30, 44 and 70 GHz) the new matched matrix filters clearly 
outperform the standard matched filters for all these quality 
indicators. The matched matrix filters decrease significantly 
the noise level, what translates into a lower detection thresh- 
old and a reduced number of false detections. The flux es- 
timation is consequently improved, with a lower dispersion 
around the true input value and a lower flux at which the 
well-know selection Eddington bias occur. The improvement 
is particularly evident at 44 GHz. At 100 GHz, however, the 
performance of the two filters is very similar. We indicate 
some possible reasons for this behavior, based on general 
analytical considerations about the structure of the filters. 

One of the most interesting ways to compare two cata- 
logs obtained with different methods is to set a fixed level of 
reliability, cut the catalogs at the corresponding points and 
compare how many detections there are in each of them. We 
find a noticeable increment of the number of true detections 
for a fixed reliability level obtained with the matched matrix 
filters with respect to the standard matched filters. In partic- 
ular, for a 95% reliability we practically double the number 
of detections at 30, 44 and 70 GHz: the ratio between the 
number of detections obtainted with the MTXF and the MF 
for these channels are 1.8, 2.2 and 1.7, respectively. 

We would like to stress once more the importance of in- 
cluding mult i- frequency information in this approach. The 
new matched multi-filters can be used to increase very sig- 
nificatively the number of extragalactic point source detec- 
tions in upcoming CMB experiments such as Planck as well 
as in current experiments such as WMAP. A work on the 



application of this technique to the 5 year WMAP data is 
in preparation. 

Although here we have tested the MTXF to the particu- 
lar case of the detection of radio sources in the low frequency 
channels of Planck, the same technique can be easily applied 
to other frequency bands, or to other fields of image analysis 
where pointlike objects appear in different frames (images). 
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